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ABSTRACT 

The distinction between chaotic and regular behavior of orbits in galactic models is an 
important issue and can help our understanding of galactic dynamical evolution. In this 
paper, we deal with this issue by applying the techniques of the Smaller (and Gener- 
alized) ALingment Indices, SALI (and GALI), to extensive samples of orbits obtained 
by integrating numerically the equations of motion in a barred galaxy potential. We 
estimate first the fraction of chaotic and regular orbits for the two-degree-of-freedom 
(DOF) case (where the galaxy extends only in the (x, y)-space) and show that it is a 
non-monotonic function of the energy. For the three DOF extension of this model (in 
the z-direction), we give similar estimates, both by exploring different sets of initial 
conditions and by varying the model parameters, like the mass, size and pattern speed 
of the bar. We find that regular motion is more abundant at small radial distances 
from the center of the galaxy, where the relative non-axisymmetric forcing is relatively 
weak, and at small distances from the equatorial plane, where trapping around the 
stable periodic orbits is important. We also find that the variation of the bar pattern 
speed, within a realistic range of values, does not affect much the phase space's fraction 
of regular and chaotic motions. Using different sets of initial conditions, we show that 
chaotic motion is dominant in galaxy models whose bar component is more massive, 
while models with a fatter or thicker bar present generally more regular behavior. 
Finally, we find that the fraction of orbits that are chaotic correlates strongly with 
the bar strength. 
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1 INTRODUCTION 

Exploring the nature of orbits in galaxies constitutes a very 
important issue, not only because of the evident astronom- 
ical interest in classifying the types of orbits that exist in 
such systems, but also because orbits are needed for con- 
structing self-consistent models of galaxies. In order to study 
and understand the structure and dynamics of a galaxy it 
is necessary to identify the type of dynamics characterizing 
the motion of stars (regular or chaotic) and estimate the 
percentages of each type in the phase space of the galaxy. 

From previous works that describe galaxies and their 
star motion, it is well-known that the analysis of periodic or- 
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bits, and their stability, can provide very useful information 
about galaxy structure. Stable periodic orbits are associated 
with regular motion, since they are surrounded by tori of 
quasi-periodic motion. Thus regular orbits are trapped in 
the vicinity of the parent periodic orbit. On the other hand 
unstable periodic orbits breed chaos. If we pick an orbit in 
the immediate phase-space vicinity of a chaotic one, we find 
that the two orbits will diverge exponentially with time. 
Such chaotic orbits fill up all the phase space region that 
is available to them. Nevertheless, recent results in galactic 
dynamics show that there are chaotic orbits that can sup- 
port galaxy features, even thin fe atures like the outer parts 
of bars, spirals an d rings (Kaufm ann &: Contopoulos 199a 
Patsis et aklllQQTl: [Patsis 2006; Romero-Gomez et al.ll2006 



2007 
,2008 
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iHarsoula fc KalapotharakosI l2009l : iTsoutsis et all [iooi). 

These orbits are often called "sticky" a nd their true chaoti c 
natu re takes very long to be revealed (jContopoulosI [20o3 ) . 
lAthanassoula et al.l (|2010l ) termed these orbits "confined 
chaotic orbits" because the structures they generate are well 
confined in configuration space. There are also several recent 
related results showing that strong loc al instability does not 
mean diffusion in phase space (e .g. Giordano fc Cincottal 
12004 ICincota fc Giordanj l2008l : ICachucho et al.1 I20f0l ). 
In our next paper, paper II, we will focus on this special 
category of chaotic orbits discussing their significance from 
an observational point of view. 

Ferrers' potentials (|Ferrerslll877l ) have proven very ef- 
ficient for studying the main properties of real bars. The 
detection of periodic orbits and their stability in these 
models have been studied in detail by many researchers 
(e.g . Ath anassoula ct al. f9 83; Pfonnigcr f984; Skokos et al. 
[2002aH bl: IPatsis et al.ll2002l . l20C)3al lbl). It is well-known that 
stable periodic orbits of low period possess in their neighbor- 
hood sizeable domains of quasiperiodic (or regular) motion. 
A number of questions arises, therefore, concerning these 
regions of stability: How far from the stable periodic orbit, 
can regular motion be sustained? Where are chaotic regions 
located in phase space and configuration space? What are 
the model's parameters that favor large islands of stability 
around the main stable periodic orbits? In this paper we 
plan to address such questions and to provide some tenta- 
tive answers. 

Chaos arises in many areas of dynamical systems and 
its effect is generally related to the type of no nlinear terms 
prese nt in the model's equations of motion (see lContopouloj 
I2OO2I for a review). The dynamics of stars in galactic po- 
tentials is a particular case, where nonlinear terms play an 
important role. The distinction between regular and chaotic 
motion is not trivial and becomes more subtle in systems 
of many DOF. Thus, it is necessary to use fast and precise 
methods to identify the nature of orbits in such models. 

Many methods have been developed over the years deal- 
ing with this problem. The inspection of successive intersec- 
tions of an orbit with a Poinca re surface of section (PSS) 
IILieberman &: Lichtenberj|l992l ) has been particularly use- 
ful for two DOF Hamiltonian systems. One of the most 
popular methods of chaos detection is the computation of 
the maximal "Lyapunov Ch aracteristic Exponent" (LCE) 
llOseleded 19681 : iBenettin et al.,1980a .b: Skokos.20fQV Along 
the same lines as LCE, several other methods of chaos detec- 
tion have been proposed based on the study of the behavior 
of deviation vector s: We may mention, e . g. the "Fast Lya- 
puno v Indicator" (|Froeschle et al.l 1 19971 : iFroeschle fc Legal 
Il998l ) and the "Mean Exponential Growth of Ne arby Or- 
bits" (ICincota fc Simdl2000l : ICincota et aLll2003l ). On the 
other hand, there also exist methods based on the anal- 
ysis of time series constructed by the coordi nates of each 
orbit such as the " Frequency M ap Analysis" (|Laskaij|l990l : 
iLaskar et all 1 19921 : lLaskailll99j 'l. More details about these 
and other re l ative methods can be fou nd in the boo k of 
IContopoulo3 (|2002l ) and in the review of lSkokoj (|2010l ) . 

In the present paper, we use a method called 
the "Smaller ALingment Index" (SALI), based on 
the properties of two deviation vectors of an orbit 
for the quick and efficient dis tinction of ch aotic mo- 
tion, originally introduced by ISkokosI (|200ll ). It has 
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Macek et al. 2010). In every case, it was confirmed to 
be a fast and reliable indicator of the chaotic or ordered 
nature of orbits. We also use SALI's recent generalization, 
the so-called "Generalized ALi gnment Index" (GALI) 
introduced bv lSkokos et"aLl l|2007h . using a set of p initially 
linearly independent deviation vectors of the system. 
Thus, following more deviation vectors, we manage to 
acquire more extra information about the complexity of 
the regular motion, i.e. the dimensionality of the invariant 
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Sections [5] and [3] 
the Lya- 



as 

detectors we use, i.e 
punov spectra and the SALI/GALI methods. In Section |4j 
we present in detail the mo del, which is co mposed by a bulge, 
a disk and a Ferrers bar (|Ferrers|[l877l ). In Section H we 
present our results for the two DOF restriction of the gen- 
eral model, calculating the regular and chaotic orbits and 
constructing also charts of the different types of motion in 
the associated phase space. Section [6] is dedicated to the 
study of the full three DOF model in terms of computing 
percentages of regular and chaotic orbits, selecting different 
sets of initial conditions and varying the parameters of the 
bar component. Thus, with the help of the SALI method, 
we first distinguish the true nature of the studied orbits in 
phase space (Section (Tjl and correlate the bar force with the 
presence of larger amount of chaotic motion in phase space. 
Finally, in Section [S] we summarize our conclusions. 



2 LYAPUNOV EXPONENTS 

The Lyapunov Characteristic Exponents (LCEs) or Lya- 
punov Characteristic Numbers (LCN) or more simple 
Characteristic Exponents are very important for the 
study of dynamical systems, for distinguishing between 
regular and chaotic b eh avior of orbits in phase space 
llBenettin et al.l Il980allbl: iLieberman fc Lichtenberd 1 19921 : 
|Pettinill2007l : ISkokos|[2010l '). In practice, the LCEs describe 
the rate of separation of infinitesimally close trajectories. 
The mathematical definition of LCE relies on Oseledecs mul- 
tiplicative theorem (|Oselededll968l ). 

A flow x(t) generated by an autonomous first-order sys- 
tem is given by: 



dx(t) 
dt 



(1) 



where F is its velocity field. Let us consider a trajectory in 
M-dimensional phase space together with a nearby trajec- 
tory, with initial conditions xo and xo -I- Axo, respectively. 
These evolve with time yielding the tangent vector Ax(xo, t) 
with its Euclidean norm: 



d(xo) = ||Ax(xo,t)|l 



(2) 
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Writing Ax = {Axi, Axm) = w, the time evolution for 
w if found by linearizing ([l]), to obtain the variational equa- 
tions: 



dt 



= J(x(^))^ 



(3) 



where J(x(t)) — dF/dx is the Jacobian matrix of the -F(x). 
The mean exponential rate of divergence of two initially 
close trajectories is: 

d(xo,f) 



(t(xo, w) = lim (-) In ■ 



(4) 



't' d(xo,0)' 

It can be shown that a exists and is finite. Furthermore, 
there is an M-dimensional basis ii of w such that for any 
w, a takes one of the M (possibly non-distinct) values 
cri(xo) = (T(xo,ei), Vi = 1,2, ...,M, which are the Lya- 
punov characteristic exponents. These can be ordered by 

size (Tl ^ (72... ^ (TM. 



3 THE METHOD OF THE SALI (AND GALI) 
SPECTRA 

Let us consider a Hamiltonian flow of A'^ DOF, an orbit 
in the 2A''-dimensional phase space with initial condition 
P(0) = (a;i(0), 2:2(0), a:;2]v(0)) and two deviation vectors 
wi(0), W2(0) from the initial point P(0). In order to com- 
pute the SALI for that orbit one has to follow the time evo- 
lution of the orbit itself as well as the two deviation vectors 
wi(f), W2(t) which initially point in two arbitrary directions. 
The evolution of the deviation vectors is given by the varia- 
tional equations ^ of the flow. At every time step the two 
deviation vectors wi (t) and W2 (t) are normalized by setting: 



W,;(f) 



W,(f) 



i = 1,2 



(5) 



iiw,(t)r 

and the SALI is computed as (|Skokoj|20"oil '): 

SALI(f) = mm {||wi(f) -I- W2(f)|| , \\wi{t) - W2(t)||} . (6) 

The properties of the time evolution of the SALI rapidly 
distinguish between ordered and chaotic motion as follows: 
The SALI fluctuates around a non-zero value for ordered or- 
bits, while it tends exponentially to zero for chaotic orbits 
l|Skokos et al.l f2003a bV In general, two different initial devi- 
ation vectors become tangent to different directions on the 
torus, producing different sequences of vectors, so that the 
SALI always fiuctuates around positive values. On the other 
hand, for chaotic orbits, any two initially different deviation 
vectors in time tend to align in the direction defined by the 
maximal LCE (mLCE). Hence, they either coincide with 
each other, or become opposite, which leads to the SALI 
falling exponentially to zero. Thus, this completely different 
behavior of the SALI helps us distinguish between ordered 
and chaotic motion in Hamiltonian systems of any dimen- 
sionality. An analytical study of SAL I 's beh avior for such 
orbits was carried out in ISkokos et al] l|2004l ). where it was 
shown that SALI a g-('^i-°"2)t^ ai,a2 being the two largest 
Lyapunov exponents. 

The usual technique to decide whether an orbit can be 
called chaotic or regular is to check, after some time inter- 
val, if its SALI has become less than a very small threshold 
value. In the following we will take this value to be equal 
to 10~*. Depending on the location of the orbit, this limit 



can be reached more or less fast, as there are phenomena 
that can hold off the final characterization of the orbit, and 
certain orbits behave as regular for long times before finally 
drifting away from regular regions and starting to wander in 
a chaotic domain. 

The generalized alignment index of order p (GALIp) is 
determined through the evolution of 2 ^ p ^ 27V initially 
linearly independent deviation vectors Wi(0),i = 1,2, ...,p, 
so it is related to the computation of many LCEs rather 
than just the maximal one. The evolved deviation vectors 
Wiit) are normalized every few time steps in order to avoid 
overfiow pro blems, but th e ir dire ctions are left intact. Then, 
according to lSkokos et al.l l|2007l ). GALIp is defined to be the 
volume of the p-parallelogram having as edges the p unitary 
deviation vectors 'Wi(t),i = 1, 2, ...,p: 



GALIp(f) =11 wi{t) A W2(t) A ... A Wp(f) 



(7) 



least two of the deviation vectors become linearly dependent, 
the wedge product in Eq. ((7| becomes zero and the GALIp 
vanishes. 

In the case of a chaotic orbit, all deviation vectors tend 
to become linearly dependent, aligning in the direction de- 
fined by the maximal Lyapunov exponent and GALIp tends 
exponentially to zero following the law: 

-[(<Ti-CT2) + (CTi~tT3) + ... + (cri-<Tp)]t 



GALIp (f) 



(8) 



where ai > . . . > ap are approximations of the first p largest 
Lyapunov exponents. In the case of regular motion, on the 
other hand, all deviation vectors tend to fall on the A'^- 
dimensional tangent space of the torus, where the motion is 
quasiperiodic. Thus, if we start with p ^ TV general deviation 
vectors, these will remain linearly independent on the A'"- 
dimensional tangent space of the torus, since there is no par- 
ticular reason for them to become aligned. As a consequence, 
GALIp in this case remains practically constant for p !^ N. 
On the other hand, for p > N, GALIp tends to zero, since 
some deviation vectors will eventually become linearly de- 
pendent, following power laws that depend on the dimension- 
ality of the torus. According t o Ichristodoulidi fc Bountia 
(2006^ and ISkokos et all (120081 ) one obtains the following 
formula for the GALIp, associated with quasiperiodic orbits 
lying on fc-dimensional tori (where k is potentially equal to, 
or smaller than the system's size): 



GALIp (t) 



constant, 
1 

1 



if 2 «; p < fc 
if fc < p s; 2^^ - fc 
if 2A - fc < p < 27V. 



(9) 



In the case where k = N, GALIp remains constant for 2 ^ 
p ^ A and decreases to zero as ~ l/i^^*'"^' for A < p < 2N. 
An efficient way to calculate GALIp is by multiplying the 
singular values Zi,i = 1, ...,p, computed through a Singular 
Value Decomposition procedure of the matrix formed by the 
deviation vectors Wi,i = l,...,p jAntonopoulos fc Bountia 
I2OO6 : ,Skokos et al.. .2008') : 



GALI„ 



(10) 



The method has been applied succe ssfully in several H amil- 
tonian systems like the FP U lattice |Skokos et aDl2008h and 
coupled symplectic maps (|Bountis et al.l 20091 ) for the de- 
tection not only of regular and chaotic motion but also the 
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dimensionality of the torus on which a regular trajectory lies 
on. 

Practically, the SALI is equivalent to GALI2 and the 
distinction between regular and chaotic motion in the 2D 
Ferrers barred model (two DOF) can be done with either 
one of them. For the full 3D version of the model, we also 
use GALI3, depending on the properties of the model (if its 
phase space is dominated by large chaotic or regular areas) 
and on our goals. GALI3 generally demands more CPU time, 
since it is necessary to follow the evolution of three deviation 
vectors instead of two. For the case of chaotic trajectories, 
however, GALI3 decays much faster and the calculation can 
stop well before the end of total time interval. Summarizing 
therefore, if someone had to estimate the amount of chaotic 
vs. regular regions in phase space, GALI3 would be more 
efficient for models where chaos is dominant, while SALI 
would be preferable for models with large regions of order. In 
our runs, we have used both SALI (GALI2) for the general 
description of chaotic vs. ordered regions, and GALI3 to 
follow specific regular orbits and study the dimensionality 
of the torus on which they lie. 



4 THE MODEL POTENTIAL 

The motion of a test particle in a 3D rotating model of a 
barred galaxy is governed by the Hamiltonian: 

H = ]^{pl+pl+pi) + V{x,y,z) - Qbixpy - ypx). (11) 

The bar rotates around its 2;-axis (short axis), while the x- 
direction is along the major axis and the y along the interme- 
diate axis of the bar. The Px,Py and Pz are the canonically 
conjugate momenta, V is the potential, ilb represents the 
pattern speed of the bar and H is the total energy of the 
orbit in the rotating frame of reference (Jacobi constant). 
The corresponding equations of motion are: 



X =Px 



9V ^ 

Px = + ^hPy, 



y ^Py - ^IbX, 



dv 

\lbPx 

dy 



(12) 



P^ 



dz 



The potential V of our model consists of three components: 

(i) A disc, represente d by a Miyamoto-Nagai disc 
ijMivamoto fc NagailflQTSh : 



Vd = 



GMr 



(13) 



where Mo is the total mass of the disc, A and B are its hor- 
izontal and vertical scale-lengths, and G is the gravitational 
constant. 

(ii) A bulge, which is modeled by a Plummer sphere 
l|Plummejll9m whose potential is: 



Vs = 



GMs 



■\/a;2 + + 2^ + 



(14) 



where Es is the scale-length of the bulge and Ms is its total 
mass. 



(iii) A triaxial Ferrers bar jFerrerdl 18771 ). the density p{x) 
of which is: 



p{x) = 



Pc(l - m^f ,m < 1 







, m 1 



(15) 



where pc — ^^^^ is the central density, Mb is the total 
mass of the bar and 

222 

m^ = ^--f|--H^-, a>b>c>Q, (16) 
c-^ 

with a, b and c being the semi-axes. The corresponding po- 
tential is: 

du 



Vb = — vrGafec- 



2/ NNI + l 



n + 1 A{u 



.(l-m^(n)) 



where 



2, , x^ y^ 

(") = — I- TT', — ^ ~2~r 



A\u) = {a^ + u){b^ + u){c^ +u), 



(17) 



(18) 



(19) 



n is a positive integer (with n = 2 for our model) and A is 
the unique positive solution of: 



m^(A) 



(20) 



outside of the bar (m ^ 1), and A = insid e the bar. 
The c orresponding forces are given analytically bv lPfennigerl 
l|l984h . 

This model has b een u sed extensively for orbital stud- 
ies by iPfennigeil (|l984D, [Skokos ct al. (2002ajj) and by 
iPatsis et all (|2002l . l2003al lbl). Our so-called standard (S) 
model has the following values of parameters G=l, 
f76=0.054 (54 km ■ sec'^ ■ kpc-^), a=&, b =1.5, c=0.6, A=?,, 
B=l, ts^QA, Mb=0.1, Ms=0.08, Md=0.82, both for its 
two DOF and three DOF versions. The units used are: 1 



kpc (length), 1000 km 



(velocity), 1 Myr (time). 



2 X 10" Mq ( mass). The total mass G(Ms + Mo + Mb) 
is set equal to 1. 

Our study is mainly focused on understanding the ef- 
fect of the bar on the dynamics of the model's phase space, 
by varying the bar mass Mb , its semiaxis lengths b, c and 
its pattern speed fit,. In order to do this we have con- 
sidered three more models: models C, B and M whose 
short z-semiaxis (c-parameter), intermediate j/-semiaxis 
(fo-parameter) and mass of the bar (M^-parameter) are 
twice those of the initial reference model. For each of these 
four {S, C, B and M) models we launch three different 
sets of initial conditions (distributions I , II , III). Hence, 
we consider the three standard models {IS, IIS, HIS) 
and their variations {IC, IIC, IIIC), {IB, IIB, IIIB) and 
{IM, IIM, HIM) , depending on the varying parameter. 
These sets of initial conditions will be discussed in detail 
in Section [T] The maximal time of integration of the orbits 
through the equations of motion and the variational equa- 
tions is set to be T=10,000A4"j/r (10 billion yrs), that corre- 
sponds to a time less than, but of the order of one Hubble 
time. We chose to take a maximum integration time that 
is the same for all orbits, independent of their energy. As 
a result, some of the orbits - and particularly the ones in 
the innermost regions of the galaxy - may have been cal- 
culated for a larger multiple of their own dynamical times 
than others. This choice, however, is necessary because it 
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allows comparison with, or applications to a galaxy or a 
simulation, where all information refers to one time only, 
the same for all orbits. It was therefore also adopted in 
previous orbital structure studies in galactic mode ls (e.g. 
lEl-Zant fc Shlosmanl |2002| . l2003l : IValluri et al] I2OI0I ) . This 
choice is essential in the case of some sticky orbits which may 
appear regular if calculated till e.g. 10,000Afj/r, but show 
clear signs of chaoticity when calculated for much longer 
times, i.e. orbits whose diffusion time scale is equal to many 
Hubble times. Such times, however, are irrelevant for our 
applications, and thus these orbits should not be counted as 
chaotic, even if in the long run they may become so. These 
issues will be discussed further in Paper II. 

The relative strength of the non-axisymmetric forces 
can be estimated by the quantity Qt, defined as: 



(21) 



The maximum of Qt over all radii shorter than the bar ex- 
tent is often refer red to as Qb. an d used a s a measure of the 
bar strength (e.g . Buta et al] 200 3. 2004; Laurikainen et al] 
I2OO4I : iButa et al.ll2005l : iDurbala et al.ll2009l ). From now on. 
we will for simplicity and conciseness, often refer to the 
"strong non-axisymmetric forcings" simply as "strong bars" . 
In Fig. [T]we show the above quantity for the 4 models we are 
going to study and discuss in this paper. From the quantity 
Qt it becomes clear that the increase of the mass of the bar 
component (model M : dotted line) gives rise to a "stronger 
bar" compared to the initial standard model (model S: solid 
line). On the other hand, an increase of the short 2;-semiaxis 
(c-parameter - model C: dashed line) or intermediate y- 
semiaxis (6-parameter - model B: dot-dashed line) has as 
a result "weaker bars". Hence, how do "strong bars" affect 
the relative regular and chaotic motion in the phase space of 
the full three DOF model? Before answering this question, 
let us first study and comprehend the dynamics of the two 
DOF model. 



5 THE TWO DOF MODEL POTENTIAL 

The two DOF Ferrers model is described by the Hamiltonian 
Eq. ([TT1) setting (z.v.) = (0.0). i.e.: 



Qbixpy — ypx 



(22) 



In the two bottom inset figures of Fig. [2^, we present two 
typical orbits in the {x, i/)-plane of the standard model. 
One is regular (bottom left inset), with initial condition: 
ix,y,pa:,py) = (0,-0.625,-0.201,-0.06) (orbit R) and 
the other chaotic (bottom right inset), with initial condi- 
tion: {x,y,p^,py) = (0,-0.625,-0.002,-0.24) (orbit C), 
both for the Hamiltonian value H = —0.360. Their qual- 
itatively different behavior is also shown in the PSS in 
Fig. [2^ ((y,pj,) -plane). The chaotic orbit C (gray points) 
tends to fill with scattered points the chaotic region, while 
the ordered orbit R (black points) creates a set of points 
that form a closed invariant curve, on the left part of 
the picture. In the top inset figure of Fig. [2^, we show 
the different morphologies of three regular orbits around 
the periodic orbits in the centers of the three largest is- 
lands of regular motion present in the associated PSS. The 
model's basic barred shape, in the {x, j/)-plane is mainly 
provided by trajectories around the main family of periodic 
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Figure 1. The quantity Qt as a function of radius. Its maximum 
value corresponds to Qb, which gives an estimate of the relative 
strength of the non-axisymmetric forces. We see that the increase 
of the mass of the bar component (model M: dotted line) corre- 
sponds to "stronger bar" compared to the initial standard model 
(model S: solid line). On the other hand, when we increase the 
short 2-semiaxis (c— parameter - model C: dashed line), or the in- 
termediate 3/-semiaxis (b— parameter model B: dash-dotted line) 
the bar becomes weaker. 



or bits, the so-called xl family (following the nomenclature 
of IContopo ulos & Papavannopoulos 1980) which are orbits 
elongated along the bar, i.e along the x-axis in our case, and 
inside corotation. The stability island on the right gives rise 
to elliptical-like orbits, around the x2 family, elongated per- 
pendicular to the bar, while the island on the left contains 
orbits which are elliptical-like but retrograde, very slightly 
elongated perpendicular to the bar and around the family 
x4. 

We then apply the SALI method to the orbits R and C 
and present their different typical evolution behaviors: For 
the chaotic orbit C (gray curve in Fig. [Ua), the SALI tends 
to zero (~ 10~^^) exponentially after some transient time, 
while for the regular orbit R, it fluctuates around a positive 
number (black curve in Fig. [^Jj) . The corresponding mLCE 
and (Ji, for these orbits is shown in Fig. [2};, where the (t\ for 
the regular orbit tends to zero, while for the chaotic orbit it 
tends to a positive number after a long integration time. 

This comparison shows clearly an advantage of SALI 
(GALI2), namely its ability to detect the chaotic charac- 
ter of an orbit faster than by simply following the mLCE, 
which typically takes a long time to converge. Even in our 
two DOF model, SALI starts to decrease exponentially af- 
ter a relatively short time (in Fig. [2]d, t ^ 10'^) while the 
corresponding mLCE cri starts to converge to some positive 
non-zero value for t ^ 10'^ (see gray curve in Fig. [2};). 

Exploiting the efficiency of SALI, we take initial con- 
ditions on the {y,Py)-\Aane (with a; = 0) and calculate the 
values of the index to detect very small regions of stability 
(or instability) more globally. We are thus able to construct 
a map of chaotic and regular regions, very similar to what 
is depicted in a PSS, but with more accuracy and higher 
resolution (albeit at a somewhat longer CPU time). Fur- 
thermore, as a byproduct of our application, we obtain an 
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Figure 2. Analysis of orbits in the two DOF case: a) Poincare surface of section in the (y, p^)— plane for a regular (R-black color) and a 
chaotic (C— gray color) orbit {H = —0.36). The points of the successive intersections of the regular orbit R with this plane create a closed 
curve, while the points of the chaotic orbit C fill with scattered points all its available region of motion. The orbits can be seen in the 
two inset figures in the bottom of the panel. In the top inset figure, we show examples of the three different morphologies of the regular 
orbits around the three main periodic orbits in the centers of the three main islands of stability. Note that the basic barred shape in the 
(i, 3/)— plane is mainly provided by trajectories from the central stable region, around the family of periodic orbits xl (orbits elongated 
along the bar's major axis). The other two stability islands give either elliptical-like orbits elongated perpendicular to the bar orbits (right 
island - family x2), or orbits which are elliptical-like but retrograde with respect the bar's pattern speed (left island - family x4). b) The 
corresponding evolution of the SALI for the R and C orbits of panel a. For the chaotic one C (gray line), the SALI decays exponentially 
fast to zero, while for the ordered one R (black line) it fluctuates around a non-zero number, c) The corresponding evolution of the 
maximal Lyapunov exponents a\ for the orbits R and C. For the chaotic orbit C (gray line) this tends after some transient time to a 
non-zero value, while for the regular orbit R (black color) it tends linearly to zero. Note that the Log(SALI) in j/-axis of panel b ranges 
from -16 to 2, while the Log((7i) in panel c ranges from -5 to 0. 



accurate estimate of the percentage of regular to chaotic or- 
bits on a surface of section of the given energy. 

For example, choosing different values of the energy and 
using a sample of 50,000 initial conditions equally spaced 
on the same (t/,py) -plane, we plot on the left column of 
Fig.[3]the PSS for H = -0.360, -0.335, -0.300, -0.260 and 
on the right the corresponding final SALI values obtained 
from the selected grid of initial conditions, where each point 
is colored according to its SALI values at the end of the 
integration. In the SALI plots, the fight gray color cor- 
responds to regular orbits, the black color represents the 
chaotic orbits/regions, while the intermediate gray shades 
between the two represent orbits with small rate of local 
exponential divergence, and/or orbits whose chaotic nature 
is revealed after long times, like for example the so-called 
sticky orbits, i.e. orbits that "stick" onto quasiperiodic tori 
for long times. Note that the fraction of these orbits (which 
lie mainly around the borders of the islands of stability) is 
very small (few % of the total amount of initial conditions) 
and can be discerned by eye in Fig. |3] only if one focuses on 
these particular regions. Thus, calculating from all tfiese ini- 
tial conditions the percentage of regular orbits, we are able 
to follow how this fraction varies as a function of the total 
energ y. This problem was alread y addressed for four mod- 
els by I Athanassoula et al.l (|l983l ) but with fewer orbits, i.e. 
larger error bars. In the present model, the chosen energy 
values start from a regime of total order, sX H = —0.46, up 
to values near the escape energy, Hesc ~ —0.20. The distinc- 
tion between chaotic and regular orbit is that SALI< 10~* 
for chaotic orbits and ^ 10~* for regular ones. In the latter 
range one does, of course, include "sticky" chaotic orbits as 
well. As is clear from Fig. |4l although the percentages of the 
regular orbits decreases sharply as the energy grows above 



100% 100% 




-0.48 -0.44 -0.40 -0.36 -0.32 -0.28 -0.24 -0.20 

Energy 

Figure 4. Percentages of regular orbits for several values of 
the energy {H = -0.460,-0.435,-0.410,-0.360,-0.335,-0.300,-0.260,- 
0.240) in the two DOF Ferrers model. Note that while the fraction 
of the regular orbits decreases sharply as the energy grows above 
H = -0.4, this tendency is reversed at higher energy values. 



H = —0.4, this tendency is reversed at higher energy values. 
Clearly, for H > —0.3, despite the fact that many stability 
islands have vanished, the size of the ordered domain in- 
creases because the main big island on the left of the PSS 
(see Fig. [3f,h) which corresponds to retrograde orbits, has 
become larger. 
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Figure 3. The PSS (left column) and SALI (right column) methods are in good agreement when surveying the same projection of the 
phase space. In the left column we show the PSS for the two DOF model potential with H = —0.360 (panel a), H = —0.335 (panel c), 
H = —0.300 (panel e) and H = —0.260 (panel g). In the right column we present regions of different values of the SALI for 50,000 initial 
conditions on the (?/, py)— plane for the same values of the Hamiltonian (panels b,d,f and h, respectively). The light gray colored areas 
correspond to regular orbits, while the dark black ones to chaotic. Note the excellent agreement between the two methods as far as the 
gross features are concerned, as well as the fact that the SALI can easily pick out small regions of stability which the PSS has difficulties 
detecting, like those around the central islands in panels b,d,f and h. 
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6 THE THREE DOF MODEL POTENTIAL 

We now turn to the three DOF model and begin a study 
of ordered and chaotic domains. Before discussing the more 
general results about the study of the phase space of the 
model, we first present in detail the behavior of our chaos 
indicators for some typical orbits. Let us choose a regular 
and a chaotic orbit and compare the efficiency of the GALI 
method with that of the mLCE's. 

In Fig. [S^,b we show the behavior of a chaotic 
orbit (Cf), with initial condition: {x,y,z,px,Py,Pz) = 
(0.5875,0.0,0.33333,0.0,0.2,0.0) and of a regular or- 
bit (Rl), with initial condition: {x,y, z,Px,Py,Pz) = 
(0.97917, 0, 0.04167, 0, -0.17778, 0). For the chaotic one, the 
GALI2 (or SALT) has become almost zero fort~ 10''. How- 
ever, as it is clear in Fig. [S^, the higher order GALI^, k = 
5, 6, indicate the chaoticity of the orbit already by t ~ 10"^, 
since they have reached very small values at that time. Note 
the good agreement between the predicted slopes related 
to its two largest Lyapunov exponents (cti ~ 0.00910 and 
(72 ~ 0.00345) given by Eq. (©. In Fig.[S]3 we show the GALI 
indices for the regular orbit Rl, where both GALl2,3 oc con- 
stant and the GALl4,5,6 decay following the power laws: 

GALl4(t) oc \, GALl5(t) oc \, GALl6(t) oc \, (23) 

t° 

obtained from Eq. ([9} for m = 3, indicating regular motion 
which lies on a 3D torus. 

One example of a low dimensional motion (lower than 3) 
is provided by the regular orbit (R2), with initial condition: 
{x,y,z,px,Py,Pz) = (0.5875,0.0,0.29770,0.0,0.33750,0.0). 
In Fig. [5j;, we show the behavior of the Log(GALI) indices 
and their slopes for the regular orbit R2. Note that only 
GALI2 remains constant while the GALl3,4,5,6 tend to zero 
following power laws predicted by Eq. Q: 

1 



GALl2(t) oc const.. 
1 



GALl5(t) oc 



GALl3(t) oc 
GALl6(f) oc 



GALUt) oc , 



(24) 



for m = 2. This implies that the orbit's motion lies on a 2D 
torus even though in three DOF Hamiltonian systems one 
generally expects the dimension of the torus to be three. 

In Fig. [6] we show the corresponding {x,y) and {x,z) 
projections for the above CI (1st column), Rl and R2 orbits 
(2nd and 3rd column respectively). The chaotic one fills up 
its available regime. Regarding the two regular ones, there 
is a clear difference reflected in their projections. The "com- 
plexity" of Rl (quasiperiodic motion on a 3D torus) is more 
pronounced than the one of the orbit R2 (quasiperiodic mo- 
tion on a 2D torus). Thus, in such cases GALI3 offers an 
extra advantage in helping us detect different "degrees" of 
regularity and does not serve only to distinguish between 
chaotic and regular motion. 
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Figure 5. Distinguishing chaotic from regular motion with GALI 
method, a) Exponential decay of all GALIs for the chaotic orbit 
CI (lin-log scale) with the predicted slopes related to its two 
largest Lyapunov exponents (cri 0.00910 and (T2 ~ 0.00345). b) 
Slopes of GALIs for the regular orbit Rl showing that its motion 
lies on a 3D torus, where GALl2,3 oc constant while the GALl4_5 g 
decay following a power law. c) Slopes of GALIs for the regular 
orbit R2 lying on a 2D torus with only GALI2 oc constant. 



7 THE DISTRIBUTION OF ORBITS IN 
PHASE SPACE 

In this section we focus on the detection and quantification 
of the different kinds of orbital motion (regular and chaotic) 
in phase (and configuration) space, as some parameters of 
the bar componen t vary. A similar study , in different poten- 
tials, was done in lELZant fc Shlosmanl l|2002l . |2003| ) , using 



the Lyapunov spectra as a tool. Here we concentrate mainly 
in the bar component and in understanding how its differ- 
ent properties can affect the general stability of the system. 
We use a relatively large number of trajectories (50,000), 
more than one set of initial conditions for the survey of the 
phase space, and the SALI method as a chaos detector for 
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Figure 6. Two projections {{x,y) and (x, z)— planes) of the ciiaotic orbit CI (1st column) and the quasiperiodic orbits Rl (2nd column) 
and R2 (3rd column). Note the different "complexity" between Rl, whose regular motion lies on a 3D torus, and R2, which lies on a 2D 
as GALI indices detected. 



the reasons discussed in the previous sections. We discuss 
the relative fraction of regular and chaotic motion not only 
as a function of its spatial location but also as a function of 
the total energy and to correlate it with the bar strength, 
i.e. the relative non-a:xisymmetric forcing. 



7.1 Initial conditions 

We now need to define the sample of orbits whose properties 
(chaos or regularity) we will examine. The best of course 
would have been to draw these from a distribution function 
of the model. This, however, is not available, so we choose 
initial conditions on different grids in phase space or from 
the density distribution and the available energy range in 
the model and we measure the corresponding fraction of 
different kinds of motion. These percentages can depend on 
the choice of the sets of initial conditions and a priori are 
not expected to be equal. For this reason, before varying any 
other parameter we should consider as carefully as possible 
what initial conditions to choose and how well this choice 
"spans" the allowed phase system of the system. 

We first considered sets of orbits in such a way as to 
favor the dynamics near the main family of periodic orbits 
xl (building blocks of the bar). To this end, we launched 
initial conditions along the bar's major x-axis with posi- 
tive momenta py and position (z) or momentum values (pz) 
in the ^-direction (distributions / and // below). We also 
considered initial conditions which follow the model's total 



Table 1. Varying parameters for all sets of initial conditions. 



distribution /model 


Mb 


b 


c 


IS, IIS, HIS 


0.1 


1.5 


0.6 


IC, lie, I lie 


0.1 


1.5 


1.2 


IB, I IB, IIIB 


0.1 


3.0 


0.6 


IM, IIM, IIIM 


0.2 


1.5 


0.6 



density p in positions (x, y, z) while their positions are arbi- 
trary (distribution /// below). More specifically, the three 
different classes (distributions) of initial conditions we use 
here are: 

• distribution /: 50,000 orbits equally spaced in the space 
{x,z,py)wit\ix e [0.0, 7.0], [0.0, 1.5], Pa e [0.0,0.45] and 

(y,p.,p.) = (0,0,0). 

• distribution II: 50,000 orbits equally spaced in the 
space {x,py,pz) with x G [0.0,7.0], Py G [0.0,0.35], p^ G 
[0.0,0.35] and {y,z,p^) = (0,0,0). 

• distribution ///: 50,000 orbits whose spatial coordi- 
nates are chosen randomly from the den sity distribution o f 
our model, using the rejection method (|Press et al.l 1 19861 ) 
within the rectangular box —a ^ x ^ a, —b ^ j/ 6, 
— c ^ z ^ c. We then draw a random value of the Hamilto- 
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Figure 7. Percentages of regular orbits as a function of (a) the energy, (b) the averaged < |2| > value, (c) the initial spherical radius 
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I, and (d) the averaged < R 
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> radius of the orbit's evolution, for models IS, IB, IC and IM. 



nian H (total energy) in the range [0, Hmax], where Hmax 
is 1.1 times the escape energy Hesc = — 0. 2CE] defined by the 
zero velocity curve going through the L\ or L2 Lagrangian 
point. Then we set {py,Pz) = (0, 0) and we calculate the px~ 
momenta by the relation: px = H{x,y, z,py,pz) (keeping the 
positive solution). Note that in this case instead of giving 
Pj,-momenta we tried px- This choice allows us to check how 
a different way of populating the phase space may affect the 
relative fraction of regular and chaotic motion. This third 
distributions has a better representation of the density, but 
still is arbitrary with respect to the velocities. 

Details of our distributions/models can be found in Ta- 
ble 1. In the next sections we start by studying the dis- 
tribution of regular and chaotic motion in phase space, for 
different choices of initial conditions, different bar parame- 
ters (sizes, mass and pattern speed). We also briefly study 
the distribution as a function of energy and of location in 
configuration space. 



^ The reason we choose an extended range of energies is because 
we wish to study some orbits that are in this range but don't 
necessarily escape during our adopted integration time. 



7.2 Percentages of regular orbits vs. energy 

We next examine the variation of the percentages of regular 
and chaotic orbits as a function of the energy and we show 
the results in Fig. [7^. We first divided the energy interval 
into 30 subintervals, each containing an equal number of or- 
bits. In every subinterval we calculated the percentages of 
regular and chaotic orbits. In Fig. [7^ we plot the percentage 
of regular orbits for energies up to ~ —0.15, i.e. up to 
a value somewhat beyond the escape energy for all models 
under study. Generally, we observe that the fraction of reg- 
ular orbits decreases as the energy increases. An important 
peak, however, occurs before the escape energy, so that for 
{H > —0.20) the fraction of regular orbits increases again! 
This non-monotonic behavior is related to the growth of 
the islands around some basic stable periodic orbits in the 
phase space and is similar to what was observed in the two 
DOF model, comparing the fraction of regular orbits for the 
common energy interval, i.e. up to Ji" ^ —0.24 (see Fig. |4]). 

Picking another set of initial conditions we find similar 
qualitative results with some small quantitative differences. 
Hence, we find the same behavior when studying the models 
of distribution //, i.e. regular motion is dominant for small 
generally energies, then a decrease and then a peak. Dis- 
tribution III manifests the same trend for small energies 
while the peak is less intense compared to the other two 
distributions. 
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Figure 8. Slices of the (x,py)-plane for different z values. Dark regions correspond to chaotic initial conditions, light gray to regular 
ones and intermediate colors to orbits with a small rate of local exponential divergence, and/or orbits whose chaotic nature is revealed 
after long times, like for example the so-called sticky orbits. The gray-scale bar represents the values of the SALI in logarithmic scale. 



7.3 Spread of regular orbits in configuration space 

We also explored the way in which regular and chaotic orbits 
are distributed along the ^-direction of the configuration 
space. Following the evolution of each orbit, we calculated 
the mean absolute value of their z-coordinate (< l^j >) and 
plot in Fig. [JJj the fraction of regular orbits as a function of 
the < \z\ > for the models of distribution I. It clearly fol- 
lows from these results that the intervals relatively near the 
{x, t/)-plane (< |z| > < 0.35) contain mainly regular orbits, 
while those with larger values of mean absolute distance z 
are mostly chaotic. Similar qualitative results are obtained 
for distribution // (not shown here). Checking distribution 
III we find good match for the relative fraction of regular 
motion for relatively small < \z\ > (up to 0.4) with the other 
models. For intermediate values (0.4 < < jzj > < 1.5) the 
motion is mostly chaotic while we don't have orbits reaching 
larger values, which is due to the way the initial conditions 
are given. 

We also looked at these percentages as a function of 



the initial spherical radius {Rsphericai = y^x^+y^ + Zq) 
and the mean spherical radius over the evolution of the 
orbits (< Rsphericai >)• Again, dividing the total range 
of the Rsphericai values in 30 subintervals, we calculated 
the percentages of regular orbits as a function of the mean 
Rsphericai valuc iu cach subintcrval. We find that the frac- 
tion of regular orbits for all models decreases sharply with 
increasing Rsphericai up to Rsphericai < 1-5, whcrc it reaches 
a minimum (see Fig.[71;). For 1.5 < Rsphericai < 7 it starts to 
increase gradually almost monotonically for model 75* and 
with some rather week fluctuations for IC. Note that the 
model IB possesses a relatively large fraction of regular mo- 
tion within the interval 2.5 < Rsphericai < 4.5 as well. Model 
lAI follows a similar trend, but with a different scale than 
IB, containing its larger fraction of regular motion in the 
same roughly intervals, while for larger Rsphericai the mo- 
tion is mainly chaotic. This result is in good agreement with 
the results in Fig. [7}1, where the horizontal axis corresponds 
to the value of the mean spherical radius over the evolution 
time of the orbits. As previously, distribution II confirms 
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well this behavior while distribution III is in good agree- 
ment for relatively small radial distances from the center of 
the galaxy and diversifies for larger ones. 

7.4 Fraction of regular motion as a function of 
distance from the equatorial plane 

We then examine whether the orbital motion is more "reg- 
ular" near the (x, y)-plane or far from it and how this re- 
lates to their behavior as a function of < \z\ > discussed 
in the subsection 17.31 and Fig. [TId. For this, we took the set 
of initial conditions from distribution / (version IC in this 
example) with initial conditions given in the space (a:;, z,py) 
and create a mesh in the (a;, p^) -plane with different slices 
in the z-direction that start from z = up to 2 = 1.25 with 
step = 0.25. In Fig. [S] we see that as the distance of the 
initial conditions from the equatorial plane increases (from 
top to bottom row and from left to right panels), large is- 
lands of stability in the (a;,pt,)-plane start to shrink, when 
X is between and 2. This trend is almost monotonic up to 
the z-slice z ~ 1 (Fig.[Sl third row, left panel). For z = 1.25 
(Fig. [51 third row, right panel) the stability region, lying be- 
tween X — 1.5 and x = 2 starts to grow. This result turns 
out to be in good agreement with the one in Fig. [TId, where 
we plot the percentages of regular orbits as a function of 
< \z\ >. This trend, as expected, is similar for the other 
two different sets of initial conditions (distributions IIC and 
I lie), since the model parameters remain the same. 

7.5 Percentages of regular orbits vs pattern speed 

Furthermore, choosing one class of initial conditions (distri- 
bution IS), we compare different models for several values 
of the bar's pattern speed fib, keeping the remaining param- 
eters constant. In this study, we try to explore the effect of 
the value of the pattern speed on the percentages of globally 
ordered or chaotic behavior of the system. 

Based on the s tructu re and crowding of the periodic or- 
bits, IContopoulo^ jlQSd ) showed that bars have to end be- 
fore corotation, i.e. that vl > a, where r_L the Lagrangian, 
or corotation, radius. A similar result was reached from the 
calculation of the s elf-consistent response to a bar forcing 
(j Athanassoula|[l980l ) . These two results are very useful, since 
they set a lower limit to the corotation radius, which can 
thus not be smaller than the bar length, but give no in- 
formation on an upper limit. Comparing the shape of the 
observed dust lanes along the leading edges of bars to that 
of the shock loci in hydro dynamic simul a tions o f gas flow in 
barred galaxy potentials, lAthanassou li (|l992al lbl) was able 
to set both a lower and an upper limit to the corotation ra- 
dius, namely tl = (1.2 ± 0.2)a. This restricts the range of 
possible values of the pattern speed for our model, from 
Qb = 0.0367, that corresponds to the Lagrangian radius 
tl ~ 1.4o, to Qb = 0.0554, that corresponds to vl ~ 1.0a. 
Using these extreme values, and the three intermediate fre- 
quencies: fib = 0.0403, fib = 0.0444 and fib = 0.0494, that 
correspond to the Lagrangian radii tl = 1.3a, tl = 1.2a 
and vl = 1.1a, respectively, we investigated how the value 
of the pattern speed affects the system and found that as 
fib increases the percentage of the regular orbits decreases 
relatively weakly. The corresponding results are given in Ta- 
ble 2. It turns out that the variation of the pattern speed 



Table 2. Models with different pattern speeds fib- 



model 


Qb 


rL 


% of Regular orbits 


mi 


0.0367032 


1.4a 


28.39 


m2 


0.0403014 


1.3a 


27.63 


ma 


0.0444365 


1.2a 


27.26 


m4 


0.0493654 


1.1a 


26.55 


ms 


0.0554349 


1.0a 


25.55 



does not affect drastically the relative fraction of regular 
and chaotic motion. One would probably need to try more 
extreme values of fib, beyond the upper and lower limits of 
the corotation radius, in order to see a significant change. 
Such values, however, would be unrealistic for real galaxies. 



7.6 Fraction of regular and chaotic trajectories 

In this section and in Fig. [9l we investigate the amount of 
regular motion in the phase space for all the three sets of 
initial conditions (distributions /, // and ///), varying the 
axial ratios {b/a or c/a parameters) and the mass (Mb pa- 
rameter) of the bar component, as described in Table 1. In 
all distributions we used 50,000 initial conditions and we 
employed the bootstrap-method (jPress et al.lll986h on sev- 
eral subsamples to make sure that this number is sufficiently 
high to provide reliable information. Computing the distri- 
butions of regular and chaotic orbits for the above classes of 
initial conditions we find the following: 

7. 6. 1 Percentages for distribution I 

Before starting with the full 3D model, we first measured the 
percentages of regular orbits in the 2D model. We use initial 
conditions equally spaced in {x,py), with {y,Px) ~ (0,0), 
X G [0.0, 7.0] and py G [0.0, 0.45]. We find relative fractions 
of regular motion equal to 45.40% for model IM, 53.43%, 
for model IS, 92.43% for model IB and 64.43% for model 
IC. This shows that at least in the 2D case, the massive 
bar model {IM) has the highest fraction of chaotic orbits, 
followed by the standard distribution IS. The two models 
with the extended short (IC), or intermediate axis (IB) 
have the least chaos. 

Let us now turn to the full 3D orbital coverage. In 
Fig. [9^ we present the percentages of the regular orbits for 
the various sets of initial conditions of the distribution /. 
Comparing the percentages, we see that increasing the mass 
of the bar in model IM increases chaotic behavior. This con- 
firms the results obtained above and in I Athanassoula et al.l 
(|1983| ) for the two DOF case. On the other hand, when the 
bar is thicker (model IC), i.e. the length of the z-axis larger, 
the system becomes more regular. Finally, the correspond- 
ing results for a fatter bar, i.e. with larger j/-axis (model 
IB) demonstrate that the increase of the intermediate axis 
of the bar also provides the system with more ordered be- 
havior. Comparing the 2D and 3D cases presented above, 
we note that the trends we find are the same, but that the 
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Figure 9. Percentages of regular orbits (SALI^S 10~*) for distribution / (left panel), // (middle panel) and /// (right panel) of initial 
conditions for all the different model versions, where the bar mass and the length of its intermediate and minor (vertical) semi-axes are 
varied. Generally, when the bar mass increases the model tends to be more chaotic, while when the size of the bar semi-axes (b or c 
parameters) increases the fraction of regular motion grows. The small difference in the trends between models IS IC and HIS — > IIIC 
is due to the different way in which the sample of initial conditions in phase space is chosen. 



relative fractions of regular orbits are somewhat higher in 
the 2D case. We will discuss this further towards the end of 
the section. 



7.6.2 Percentages for distribution II 

Similarly, in Fig. |9}3, we present percentages of the regular 
orbits for initial conditions selected for distribution // and 
models S, M, B and C. As before, the increase of the bar 
mass (model IIM) causes more extensive chaos, while for a 
thicker bar or fatter bar (models II B and IIC), the system 
is again more regular. It is thus clear that the results of 
distributions / and II are similar, presumably dues to the 
fact that they contain initial conditions that cover the phase 
space in a similar manner, i.e. they favor the region around 
the xl family of periodic orbits. 

7.6.3 Percentages for distribution III 

Finally, we plot in Fig. ^ the percentages of regular orbits 
for several sets of initial conditions from distribution III. 
A first general observation is that the fraction of ordered 
motion is significantly smaller than is distributions I and II, 
for all models. The basic reason for this difference is related 
to the way that the momenta of the positions are given in 
this distribution of initial conditions. In particular, all the 
orbits in this case are launched with a Px momentum instead 
of a Py as in distributions / and II. As a consequence, the 
big island of stability around the main family of periodic 
orbits xl is populate with fewer orbits and thus in general 
more chaotic orbits are launched and measured. Neverthe- 
less, again the increase of the bar mass in the HIM model 
results in more chaotic behavior. As for the IB model, a 
thicker bar in the y-direction turns out to make the sys- 
tem more regular. Regarding the case of larger z-semiaxis 
(model IIIC), we notice a small decrease in the percent- 
ages of regular orbits (9.74% now, from 10.46%). It turns 
out that this particular different distribution of initial con- 
ditions doesn't reveal quite the same trend as in the previous 
cases {IS IC and IIS IIC). 
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Figure 10. The relative fraction of chaotic motion as a func- 
tion of bar strength Q(, increases (see Fig[lJ. Filled squares (■) 
correspond to the models M, S, C and B (from right to left, 
respectively) for distribution I. while filled triangles (T) are for 
distribution // and filled circles (•) for distribution 



Concluding this section dedicated to the dynamical 
study of regular and chaotic motion in the phase space 
is that their corresponding fraction is basically related to 
the way that one populates the main stable (or unsta- 
ble) periodic orbits of the system. We find that motion 
near the plane is generally stable (non chaotic) and that 
trajectories spending the largest fraction of time at large 
z are usually chaotic (with large Lyapunov exponents). 
These results are in agreement with the results obtained by 
lEl-Zant fc Shlosman (2002, 2003 ) for different models. The 
increase of the bar mass leads to more chaotic mot ion in the 
phase space as found in lAthanassoula et al.l (|l983l ) for a two 
DOF model, while in general the increase of the bar y or 
z-semiaxis leads to more regularity. 
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7.7 The effect of bar strength 

Reviewing the above percentages of regular (and chaotic) 
motion in phase space, we see that they change in the same 
way when the basic parameters and properties of the model 
are varied. In order to understand this better we examine 
next the effect of the bar strength, measuring it as dis- 
cussed in Section [4l i.e. from the relative strength of the 
non-axisymmetric forcing. In Fig. [10] we plot the relative 
fraction of chaotic motion as a function of bar strength Qf^. 
Black filled squares correspond to distribution /, filled trian- 
gles to distribution // and filled circles to distribution III. 
This figure reveals that, provided the choice of initial condi- 
tions for the orbits is the same or similar, there is a tight cor- 
relation between the two quantities, the amount of chaos in- 
creases with increasing bar strength, as intuitively expected. 
What is unexpected, however, is how tight all these correla- 
tions are, with rank correlation coefficients ~ —0.936 for all 
distributions. 

Fig. [To] also shows that, except for the bar strength, 
the distribution of initial conditions chosen is also impor- 
tant. Distributions / and II differ little and their best fit- 
ting straight lines are roughly parallel and little displaced 
from each other. This is in good agreement with our pre- 
vious statement (subsect. I7.ip that they both cover well 
the regular orbits trapped around the stable orbits of the 
xl family, which are in fact the backbones of the bar 
([Athanassoula at al.ll 198^ 1. On the other hand, distribution 
/// has a much higher fraction of chaotic orbits, as could 
also be seen from Fig. [9];, attesting that the correspond- 
ing distribution of initial conditions is more chaotic, for the 
reasons we already described. 



8 CONCLUSIONS 

In this paper, paper I, we studied the detection and distri- 
bution of regular and chaotic motion in the phase space of 
barred galaxy models. Our results referring to the dynam- 
ical properties of the 2D/3D model with a Ferrers bar can 
be summarized as follows: 

(i) We applied successfully the GALI method to distin- 
guish both qualitatively and quantitatively between regular 
and chaotic orbits. We showed its efficiency and its advan- 
tage in detecting fast and accurately the chaotic nature of a 
trajectory. 

(ii) In 2D systems, using the SALI (GALI2) method, we 
were able to identify efficiently and fast tiny regions of regu- 
lar or chaotic motion, which are not clearly visible on PSSs. 

(iii) Using GALl2,3 we detected regular motion in low di- 
mensional tori, i.e. examples of regular orbits of the three 
DOF Ferrers' model that lie on a 2D torus, while the torus' 
expected dimension is generally three. Concerning chaotic 
orbits, the GALl3,4,5,6 indices decay faster than SALI and 
are able to detect chaos at early times well before this is 
evident from the calculation of the mLCE. 

(iv) We also tried different models and different distribu- 
tions of initial conditions for the orbits to test the fraction of 
regular and chaotic motion in various cases. We found that 
regular orbits are generally dominant at relatively small ra- 
dial distances from the center of the galaxy and at small 
distances from the equatorial plane. 



(v) We tested the effect of varying the bar pattern speed 
for our standard model and found that, within a realistic 
range of values, this parameter does not affect much the 
phase space dynamics. Within these limits, we found that 
the fraction of regular motion varies only by about 3 per 
cent, in the sense that the bars at the slow limit of the 
realistic range have more regular motion than the bars in 
the fast limit. 

(vi) We varied the values of the main parameters of our 
models, such or the mass and the axial ratio of the bar com- 
ponent (in the 3D model). Chaos turns out to be dominant in 
galaxy models whose bar component is more massive, while 
models with a thicker or fatter bar present generally more 
regular behavior although the initial conditions are given 
can in general affect somewhat the relative percentages. 

(vii) We found a very strong correlation (rank correlation 
coefficient at ~0.94) between the fraction of chaotic orbits 
and the relative strength of the non-axisymmetric forcing. 
This holds for all our three initial conditions distributions 
taken individually, even though they were specifically cho- 
sen so as to represent different types of orbits. We can thus 
conclude that strong non-axisymmetric forcings (i.e. strong 
bars) are the main cause of the presence of large amount of 
chaotic motion in phase space. 

In our next paper of this series (paper II) we focus on 
the significance of the "different degrees" of chaotic mo- 
tion, like the strong and weak chaos, and their significance 
from an astronomical point of view. The fact that many 
"sticky" /chaotic orbits often persist and behave in a reg- 
ular manner for very large time intervals, before showing 
their chaoticity, makes them astronomically important for 
supporting the galaxy structure. We present there a first 
approach in filtering these weak chaotic orbits from the 
strongly chaotic ones and classify them observationally as 
effectively ordered motion. 
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